A C°-WEAK GALERKIN FINITE ELEMENT METHOD FOR THE 
BIHARMONIC EQUATION 
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Abstract. A C°-weak Galerkin (WG) method is introduced and analyzed for solving the bi- 
harmonic equation in 2D and 3D. A weak Laplacian is defined for C° functions in the new weak 
formulation. This WG finite element formulation is symmetric, positive definite and parameter free. 
Optimal order error estimates are established in both a discrete H 2 norm and the L? norm, for the 
weak Galerkin finite element solution. Numerical results are presented to confirm the theory. As a 
technical tool, a refined Scott-Zhang interpolation operator is constructed to assist the corresponding 
error estimate. This refined interpolation preserves the volume mass of order (k + 1 — d) and the 
surface mass of order (k + 2 — d) for the Pk+2 finite element functions in d-dimensional space. 
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1. Introduction. We consider the biharmonic equation of the form 

(1.1) A 2 u = /, in Q, 

(1.2) u = g, on 5f2, 

(1.3) 7T = ^ on ^> 

on 

where is a bounded polygonal or polyhedral domain in M. d for d = 2,3. For the 
biharmonic problem (1.1) with Dirichlet and Neumann boundary conditions (1.2) and 

(1.3) , the corresponding variational form is given by seeking u € H 2 (fl) satisfying 
u\an = 9 and §^|an = <f> such that 

(1.4) (Au,Av) = (/,«) Vu G H 2 (il), 

where Hq(VL) is the subspace of H 2 (fl) consisting of functions with vanishing value 
and normal derivative on dfl. 

The conforming finite element methods for the forth order problem (1.4) require 
the finite element space be a subspace of H 2 (fl). It is well known that constructing 
H 2 conforming finite elements is generally quite challenging, specially in three and 
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higher dimensional spaces. Weak Galerkin finite element method, first introduce 
in [23] (see also [22] and [16] for extensions), by design is to use nonconforming 
elements to relax the difficulty in the construction of conforming elements. Unlike the 
classical nonconforming finite element method where standard derivatives are taken 
on each element, the weak Galerkin finite element method relies on weak derives 
taken as approximate distributions for the functions in nonconforming finite clement 
spaces. In general, weak Galerkin method refers to finite element techniques for partial 
differential equations in which differential operators (e.g., gradient, divergence, curl, 
Laplacian) are approximated by weak forms as distributions. 

A weak Galerkin method for the biharmonic equation has been derived in [18] by 
using totally discontinuous functions of piecewise polynomials on general partitions 
of arbitrary shape of polygons/polyhedra. The key of the method lies in the use of a 
discrete weak Laplacian plus a stabilization that is parameter-free. In this paper, we 
will develop a new weak Galerkin method for the biharmonic equation (1.1)-(1.3) by 
redefining a weak Laplacian, denoted by A w , for C° finite clement functions. Com- 
paring with the WG method developed in [18], the C°-weak Galerkin finite element 
formulation has less number of unknowns due to the continuity requirement. On the 
other hand, due to the same continuity requirement, the C°-WG method allows only 
traditional finite element partitions (such as triangles/quadrilaterals in 2D), instead 
of arbitrary polygonal/polyhedral grids as allowed in [18]. 

A suitably-designed interpolation operator is needed for the convergence analysis 
of the C°-weak Galerkin formulation. The Scott-Zhang operator [21] turns out to 
serve the purpose well with a refinement. This paper shall introduce a refined version 
of the Scott-Zhang operator so that it preserves the volume mass up to order (k+l—d), 
and the surface mass up to order (k + 2 — d), when interpolating H 1 functions to the 
Pk+2 C°-finite element space: 



where T is any triangle (d = 2) or tetrahedron (d = 3) in the finite element, and E 
is an edge or a face-triangle of T. With the operator Qq, we can show an optimal 
order of approximation property of the C°-finite element space, under the constraints 
of weak Galerkin formulation. Consequently, we show optimal order of convergence 
in both a discrete H 2 norm and the L 2 norm, for the C° weak Galerkin finite clement 
solution. 

The biharmonic equation models a plate bending problem, which is one of the 
first applicable problems of the finite clement method, cf. [9, 2, 10, 28]. The standard 
finite element method, i.e., the conforming element, requires a C 1 function space of 
piecewise polynomials. This would lead to a high polynomial degree [2, 26, 27, 24], or 
a macro-element [10, 6, 9, 12, 20, 25], or a constraint element (where the polynomial 
degree is reduced at inter-element boundary) [3, 19, 28]. Mixed methods for the 
biharmonic equation avoid using C 1 element by reducing the fourth order equation 
to a system of two second order equations, [1, 8, 11, 14, 17]. Many other different 
nonconforming and discontinuous finite element methods have been developed for 
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solving the biharmonic equation. Morley element [13] is a well known nonconforming 
element for its simplicity. C° interior penalty methods are studied in [5, 7, 15], which 
are similar to our C°-weak Galcrkin method except there is no penalty parameter 
here. 

2. Weak Laplacian and discrete weak Laplacian. Let I? be a bounded 
polyhedral domain in R d ,d = 2,3. We use the standard definition for the Sobolev 
space H S (D) and their associated inner products (-,-) Si d, norms || • || s ,d, and semi- 
norms | • \ s ^d f° r an y s > 0. When D = ft, we shall drop the subscript D in the norm 
and in the inner product. 

Let T be a triangle or a tetrahedron with boundary dT. A weak function on 
the region T refers to a vector function v = {vQ,v n } such that vq E L 2 (T) and 
v„ • n £ H~i{dT), where n is the outward normal direction of T on its boundary. 
The first component v can be understood as the value of v on T and the second 
component v„ represents the value Vw on the boundary of T. Note that v n may not 
be necessarily related to the trace of Vwo on dT. Denote by W(T) the space of all 
weak functions on T; i.e., 

(2.1) W(T) = {v = {« ,v„} : v EL 2 (T), v n ■ n E H~i (dT)} . 

Let (-,-)t stand for the L 2 -inner product in L 2 (T), (-,-)qt be the inner product in 
L 2 (dT). For convenience, define G 2 {T) as follows 

G 2 (T) = W: y&H\T), Aip e L 2 (T)}. 

It is clear that, for any (p g G 2 (T), we have Vip € H(div , T). It follows that V<p • n e 
H~i(dT) for any ip G G 2 (T). 

Definition 2.1. The dual of L 2 (T) can be identified with itself by using the 
standard L 2 inner product as the action of linear functionals. With a similar inter- 
pretation, for any v £ W(T), the weak Laplacian of v = {v a , v„} is defined as a linear 
functional A w v in the dual space of G 2 {T) whose action on each ip e G 2 (T) is given 
by 

(2.2) (A w v, ip) T = (wo, Aip) T - (v , Vp ■ n) 9T + (v„ • n, tp) 9T , 
where n is the outward normal direction to dT. 

The Sobolev space H 2 (T) can be embedded into the space W(T) by an inclusion 
map i w : H 2 (T) W(T) defined as follows 

» w (0) = {#r,V0|ar}, 4>eH 2 (T). 

With the help of the inclusion map iw, the Sobolev space H 2 (T) can be viewed as a 
subspacc of W(T) by identifying each (f> £ H 2 (T) with iyy (</>)■ Analogously, a weak 
function v — {v a ,\r n } e W{T) is said to be in H 2 (T) if it can be identified with a 
function <p £ H 2 (T) through the above inclusion map. It is not hard to see that the 
weak Laplacian is identical with the strong Laplacian, i.e., 

A w iw(v) = Av 

for smooth functions v E H 2 (T). 
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Next, we introduce a discrete weak Laplacian operator by approximating A w 
in a polynomial subspace of the dual of G 2 (T). To this end, for any non- negative 
integer r > 0, denote by P r {T) the set of polynomials on T with degree no more than 
r. A discrete weak Laplacian operator, denoted by A w ^t, is defined as the unique 
polynomial A w ^ : tv G P r {T) that satisfies the following equation 

(2.3) (A w , r . T v, tp) T = (v , Atp) T - (v , Vtp ■ n) 9T + (v„ • n, tp) dT V<£ e P r {T). 

Recall that v„ represent the Vv on e G dT. Define v„ = (Vv ■ n)n = v n n. Obviously, 
v„ • n = v„ • n. Since the quantity of interest is not v„ but v n • n, we can replace v„ 
by v n = v n ri from now on to reduce the number of unknowns. Scalar v n represents 
Vv • n. 

3. Weak Galerkin Finite Element Methods. Let Th be a triangular (d = 2) 
or a tetrahedral (d — 3) partition of the domain £1 with mesh size h. Denote by £h 
the set of all edges or faces in Th, and let £° = £h\dQ be the set of all interior edges 
or faces. 

Since v„ = v n n with v n representing Vv ■ n, obviously, v n is dependent on n. To 
ensure v n a single values function on e 6 4, we introduce a set of normal directions 
on £ h as follows 

(3.1) T>h = {n e : n e is unit and normal to e, e e £h}. 

Then, we can define a weak Galerkin finite element space Vh for k > as follows 

(3.2) V h = {v = {v 0l v n n e } : v G V , v n \ e G P fe+ i(e),e C dT}, 
where v n can be viewed as an approximation of Vv • n e and 

(3.3) V = {veH 1 (n) ] v\ T eP k+2 (T)}. 
Denote by V® a subspace of Vh with vanishing traces; i.e., 

Vh ={v = {v , v n n e } G V h ,v \e - 0, v n \ e = 0, e C dT n 90}. 

Denote by the trace of Vh on <9£1 from the component vo . It is easy to see that 
consists of piecewise polynomials of degree k + 2. Similarly, denote by the trace of 
Vh from the component v n as piecewise polynomials of degree k + 1. Let A^^ be the 
discrete weak Laplacian operator on the finite element space Vh computed by using 

(2.3) on each element T for k > 0; i.e., 

(3.4) {A w . k v)\ T = A WykyT {v\ T ) Vv G V h . 

For simplicity of notation, from now on we shall drop the subscript k in the notation 
A W: k for the discrete weak Laplacian. We also introduce the following notation 

{A w v, A w w) h = ^2 ( A ™ w ' A ^w) T - 
Ten 

For any Uh — {u ,u n n e } and v = {v ,v n n e } in Vh, we introduce a bilinear form as 
follows 

(3.5) s(u h ,v)= ^ h^iVuo ■ n e - u n , Vw • n e - v n ) dT . 

Ten 
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The stabilizer s(u h ,v) defined above is to enforce a connection between the normal 
derivative of u along n e and its approximation u n . 

Weak Galerkin Algorithm 1. A numerical approximation for (l.l)-(l.S) 
can be obtained by seeking Uh — {uq, u n n e } G Vh satisfying uq = Qi>g and u n = 
(n • n e )Q n (j) on dil and the following equation: 

(3.6) (A K !H, A w v) h + s(u h , v) = (/, v ) V» = {v , w„n e } g Vfi, 

where Qbg and Q n <f> are the standard L 2 projections onto the trace spaces Ah and Th, 
respectively. 

Lemma 3.1. The weak Galerkin finite element scheme (3.6) has a unique solu- 
tion. 

Proof. It suffices to show that the solution of (3.6) is trivial if / = g = <j> = 0. To 
this end, assume / = g = <f> = and take v = Uh in (3.6). It follows that 

{A w u h , A w u h ) h + s{u h , u h ) = 0, 

which implies that A. w u h = on each element T and Vu • n e = u n on dT. We 
claim that Auh = holds true locally on each element T. To this end, it follows from 
A w Uh = and (2.3) that for any ip g Pk(T) we have 

(3.7) = (A w u h , <p) T = (u , A<p) T - (u , Vtp ■ n) dT + (u n n e ■ n, f) gT 

= (Au , <p)t + (u n n e ■ n - Vu ■ n, ip) gT 
= (Au ,ip) T , 

where we have used 

(3.8) u„n e • n — Vmq ■ n = ±(u„ — Vuo ■ n e ) = 

in the last equality. The identity (3.7) implies that Awo = holds true locally on 
each element T. This, together with V«o • n e = u„ on dT, shows that Uh is a smooth 
harmonic function globally on £1 The boundary condition of u = and u n = then 
implies that Uh = on O, which completes the proof. □ 

4. Projections: Definition and Approximation Properties. In this sec- 
tion, we will introduce some locally defined projection operators corresponding to the 
finite element space Vh with optimal convergent rates. 

Let Qo : H 1 ^!) — > Vq be a special Scott-Zhang interpolation operator, to be 
defined in (A. 9) in Appendix, such that for given v g i7 1 (fi), Q v g Vq and for any 

TeTh, 

(4.1) (Qqv, Aip) T - (Q v, V<p ■ n) dT = (v, A<p) T - (v, Vtp ■ n) dT , V<p e P k (T), 
and for < s < 2 

(4.2) ( h2s W u - Q^\\It) 1/2 < C^+ 3 || U || fe+3 . 

Ten 

Now we can define an interpolation operator Qh from H 2 (il) to the finite element 
space Vh such that on the element T, we have 

(4.3) QhU = {Qqu, (Q n (Vu ■ n e ))n e }, 



G 



where Q is defined in (A. 9) and Q n is the I? projection onto Pk+i{e), for each e C dT. 
In addition, let Q/j be the local L 2 projection onto Pk(T). For any ip G Pk(T) we have 

(A w Q h u, ip) T = (Q u, Alp) t - {Qqu, Vip ■ n) 9T + (Q n (Vu ■ n e )n e • n, tp) 9T 
= (u, A(p) T - (u, V<p • n) 9T + (Vu • n, ip) dT 
= (Am, <p) T = (Q h Au, <p) T , 

which implies 

(4.4) A w Q h u = Q h (Au). 

The above identity indicates that the discrete weak Laplacian of a projection of u is 
a good approximation of the Laplacian of u. 

Let T G Th be an element with e as an edge or a face triangle. It is well known 
that there exists a constant C such that for any function g £ H 1 (T) 

(4-5) \\9fe<C(h^\\gf T + h T \\Vg\\ 2 T ). 

Define a mesh-dependent semi-norm ||| • ||| in the finite element space Vh as follows 

(4.6) IHII 2 - (A w v, A w v) h + s(v, v), v G V h . 

Using (4.5), we can derive the following estimates which are useful in the convergence 
analysis for the WG-FEM (3.6). 

Lemma 4.1. Let w G H k+3 (n) and v G Vh- Then there exists a constant C such 
that the following estimates hold true. 

(4.7) Y, K Aw - ^ Aw ' ( Vw ° - w ™ ne ) • n >^l ^ c/i fc+1 |M| fc+3 |H|, 

THT h 

(4.8) Y h^KiVQow) ■ n e - Q n (Vw ■ n e ), Vw • n e - v n ) dT \ 

T£T h 

< Ch k+1 \\w\\ k+3 \\\v\\\. 



Proof. To derive (4.7), we can use the Cauchy-Schwarz inequality, (3.8), the trace 
inequality (4.5), and the definition of Qh to obtain 

Y \( Aw - QhAw, (Vw - v n n e ) ■ n) 9T \ 
TeT h 

<\J2 MlAw - QfcAwUlr j ( h^W^vo ■ n e - v n f dT j 
\TeT h / \TeT h J 

<c(j2 (W Aw - ®hAw\\ 2 T + h 2 T \\V(Aw - Q h Aw)\\ 2 T ) J HI 
\TeT h ) 

<Ch k+1 \\w\\ k+3 \\\v\\\. 
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As to (4.8), we have from the definition of Q n , the Cauchy-Schwarz inequality, the 
trace inequality (4.5), and (4.2) that 



^ /i T 1 |((VQ w) • n e - Q n (Vw ■ n e ), Vu • n e - u n ) 



Ten 



^ /i T 1 |((VQ w) • n e - Vw • n e , Vu ■ n e - u„)st 



Ter h 



dT 



- I X! /l T 1 H( V( 3ow - Vffl) • n e || aT I I 53 /ly^lVuo • n e - v n \ 
\Ten ) \TeT h / 

<c(j2 (h T 2 \\VQow - Vw\\ 2 T + \\VQ w - Vu;||i,t)) ||v||| 
\Ten / 

<Cft fc+1 ||«;|| fc+ 3||«l|. 
This completes the proof. □ 

5. An Error Equation. We first derive an equation that the projection of the 
exact solution, QhU, shall satisfy. Using (2.3), the integration by parts, and (4.4), we 
obtain 

(A w Q h u, A w v) T 

= {vo, A{AwQhu)) T + ((v n n e ) ■ n, A w Q h u) dT - (v ,V(A w Q h u) ■ n) dT 
= {Av , A w Q h u) T + (v ,V(A w Q h u) ■ n) dT - (Vv ■ n, A w Q h u) dT 

+ ((v n n e ) ■ n, A w Q h u) aT - (v ,V(A w Q h u) ■ n) 9T 
= (Avo,A w Q h u) T - ((Vuo ~ v n n e ) ■ n, A w Q h u) dT 
= {Av Q ,q h Au) T - <(V«o - v n n e ) ■ n,Q h Au) dT 
= (Am, Au ) t - ((Vv - v n n e ) ■ n, Q h Au) dT , 

which implies that 

(5.1) (Am, Av )t = {A w Q h u, A w v) T + ((V«o - v n n e ) ■ n, Q h Au) dT . 
Next, it follows from the integration by parts that 

(Am, Au )t = (A 2 u, v ) T + (Au, Vv ■ n)g T - (V(Au) • n, v )q T . 
Summing over all T and then using the identity (A 2 w, u ) = (/, v a ) we arrive at 

53 ( A «. A «o)t = (/, «b) + ^ (Au, Vv ■ n) dT 
Ten Ten 

= (/, w o) + 53 ^ Au ' ( Vw ° - u « ne ) ■ n ^ T - 

Ten 

Combining the above equation with (5.1) leads to 

(5.2) (A w Q h u, A w v) h = (/, v Q ) + 5] (Aw - Q h Au, (Vw - w„n e ) • n) dT . 

Ten 
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Define the error between the finite element approximation Uh and the projection of 
the exact solution u as 

e h := {e , e„n e } = {Q u - u , (<9„(Vu • n e ) - u n )n e }. 

Taking the difference of (5.2) and (3.6) gives the following error equation 

(5.3) (A w e h , A w v) h + s(e h , v) = (Au - Q h Au, (Vv - v n n e ) ■ n) 9T 

Ten 

+ s(Q h u,v) WGVfc . 
Observe that the definition of the stabilization term s(-, •) indicates that 

s(Q h u, v) = h^dVQau) ■ n e - Q n (Vu ■ n e ), Vu ■ n e - w„) 9T - 

6. Error Estimates. First, we derive an estimate for the error function e/j in 
the natural triple-bar norm, which can be viewed as a discrete iJ 2 -norm. 

Theorem 6.1. Let Uh € Vh be the weak Galerkin finite element solution arising 
from (3.6) with finite element functions of order k + 2 > 2. Assume that the exact 
solution of (1.1)- (1.3 ) is regular such thatu € i? +3 (Q). Then, there exists a constant 
C such that 

(6.1) l\u h -Q h u\l<Ch k+1 \\u\\ k+3 . 



Proof. By letting v = eh in the error equation (5.3), we obtain the following 
identity 

UK ||| 2 = ^ (Au - Q h Au, (Ve - e„n e ) • n) 9T 

TGT h 

+ ^2 h^dVQou ■ n e - Q n {Vu ■ n e ), Ve • n e - e n ) dT . 
Ten 

Using the estimates of Lemma 4.1, we arrive at 

IIM! 2 < ch k+1 

which implies (6.1). This completes the proof of the theorem. □ 

Next, we would like to provide an estimate for the standard L? norm of the first 
component of the error function e^. Let us consider the following dual problem 

(6.2) A 2 w = e in fi, 

(6.3) w = 0, on dQ, 

(6.4) Vw ■ n = on 90. 

The H 4 regularity assumption of the dual problem implies the existence of a constant 
C such that 

(6.5) IHU < CIMI- 
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Theorem 6.2. Let u h e V h be the weak Galerkin finite element solution arising 
from (3.6) with finite element functions of order k + 2 > 3. Assume that the exact 
solution of (1.1)-(1.3) is regular such that u G H k+3 (fl) and the dual problem (6.2)- 
(6.4) has the H 4 regularity. Then, there exists a constant C such that 

(6.6) IIQou-uoll <C7* fc+3 ||u|| fc+3 . 

Proof. Testing (6.2) by error function eo and then using the integration by parts 
gives 

||e || 2 = (AVeo) 

= ^(Aw,Ae ) T - J2 ( Aw > Ve o ' n )dT 
Ten Ten 

= ( Au,J Ac o)t - X! (Aw, (Ve - e„n e ) • n) aT . 
Ten Ten 

Using (5.1) with w in the place of u, we can rewrite the above equation as follows 

||e || 2 = (A w Q h w, A w e h ) h - (Aw - QhAu>, (Ve - e„n e ) • n) aT . 

Ten 

It now follows from the error equation (5.3) that 

(A w Q h w, A w e h ) h = ^ ( Au ~ QhAu, (VQ w - Q n (Vw ■ n e )n e ) • n} dT 

Ten 

- s(e h , Q h w) + s(Q h u, Q h w). 
Combining the two equations above gives 

(6.7) ||e || 2 = - ( Aw - ^hAw, (Ve - e„n e ) • n) 9T 

Ten 

+ ^ (Au - QhAu, (VQqw - Q n (Vw ■ n e )n e ) • n) dT 

Ten 

s(e h , Q h w) + s(Q h u, Q h w). 

Using the estimates of Lemma 4.1, we can bound two terms on the right-hand side of 
the equation above as follows 

^2 I (Aw - QhAu), (Ve - e„n e ) • n) gT \ < Ch 2 \\w\\ 4 \\\e h \\\, 
Ten 

\s(e h ,Q h w)\ < C7i 2 ||w|| 4 ||e,J. 
It follows from (3.8) and the definition of Q n and Q that 

(6.8) \\(VQ w - Q n (Vw ■ n e )n e ) • n e \\ dT = \\(VQ w - Q n (Vw ■ n e )n e ) ■ n\\ dT 

= ||V<3ow • n - Q„(Vtu • n)\\ dT < \\VQ w ■ n - Vw • n\\ dT 
+ ||Vw • n - Q n (Vw • n)||ar < C\\VQ Q w • n - Vw • n\\ dT . 
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Using (6.8) and (4.5), we have 

(An - Qh^u, {VQ w - Q n (Vw ■ n e )n e ) • n) dT 

1/2 / \ 1/2 



or 



< C ( h W Au - ®h&u\\lj) ( Y h^WiVQow - Vw) • n| 

\T£T h ) \TeT h 

< C ( (W Au ~ ^h^uf T + 4||V(A« - Q h Au)\\ 2 T ) J • 

\TeTh / 

( E (^ 2 ll v Qo^ - VHII + ||v(vq„«> - v«;)|fr) ) 

\Ter h ) 
<C^+ 3 || M || fe+3 |h|| 4 . 

Using (6.8) and (4.5), we have 

s(QhU,Qhw) = Yj /i _1 (V(<3o«) • n e - Q n (Vu ■ n e ), V(Q w) • n e - Q„(Vw • n e ))dT 

T£T h 

< ( E ^ _1 HVQo« - Vu) • n||| r > ] C £ ft-iVQo^ - V«d) • n||^ 
\Ter h / \Ter h / 

<^ fe+3 || U |U +3 |h||4. 

Substituting all above estimates into (6.7) and using (6.1) give 

IMI 2 <^ fe+3 |MU +3 |M| 4 . 

Combining the above estimate with (6.5), we obtain the desired L 2 error estimate 
(6.6). □ 

7. Numerical Experiments. This section shall report some numerical results 
for the C° weak Galerkin finite element methods for the following biharmonic equa- 
tion: 

(7.1) A 2 u = f in ft, 

(7.2) u = g on <9Q, 

(7.3) — = v on an. 

on 

For simplicity, all the numerical experiments are conducted by using k = or k = 1 
in the finite element space Vh in (3.2). 

If <f> € PoC^O (i-e. k = 0), the above equation can be simplified as 
(A^w, 0) T = (u„n e • n, <j)) 9T . 



The error for the C°-WG solution will be measured in four norms defined as 
follows: 
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H 1 semi-norm: 



\v~voh 



EX 



Vu - \7v \ 2 dx. 



Ten 



Discrete H 2 norm: 



IIMII 2 = 51 \\A w v\\t+ h 1 W T7v • n e ~ v n\\dT, 

Ten Ten 

Element-based L 2 norm: 

WQov - v a \\ 2 = ^ / \Q v - v \ 2 dx, 
Ten Jt 

Edge-based L 2 norm: 

\\Q„(Vv ■ n e ) -v n \\l = ^2 h / \Q n (Vv ■ n e ) - v n \ 2 ds. 



ee£ h 



7.1. Example 1. Consider the biharmonic problem (7.1)-(7.3) in the square 
domain f2 = (0, l) 2 . Set the exact solution by 

U = a; 2 (l-.T) 2 y 2 (l-y) 2 . 



Table 7.1 

Example 1. Convergence rate for element PziT) — Pi(e) (k = 0). 



h 


\\u - Uo||l 


\\\u h - Qhu\\\ 


||«o - Qou\\ 


||<9„(Vu • n e ) - u n \\ b 


1/4 


6.8858c-03 


6.0250e-02 


1.4563e-03 


4.3364e-03 


1/8 


1.7465c-03 


3.0867c-02 


3.8153c-04 


1.4617e-03 


1/16 


4.3885c-04 


1.5555c-02 


9.6991c-05 


4.0941e-04 


1/32 


1.0982e-04 


7.7916c-03 


2.4350e-05 


1.0558e-04 


1/64 


2.7458c-05 


3.8972c-03 


6.0931e-06 


2.6601c-05 


1/128 


6.8645e-06 


1.9487c-03 


1.5236e-06 


6.6629e-06 


Conv.Rate 


1.9949 


9.9160c-01 


1.9829 


1.8865 



It is easy to check that 



„ du 
u\dn = 0, — = 0. 
on 



The function / is given according to the equation (7.1). 

The test is performed by using uniform triangular mesh. The mesh is constructed 
as follows: 1) partition the domain into n x n sub-rectangles; 2) divide each square 
element into two triangles by the diagonal line with a negative slope. The mesh size is 
denoted by h = 1/n. Table 7.1 shows the convergence rate for C°-WG solutions based 
on k — in four norms respectively. The numerical results indicate that the WG 
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solution is convergent with rate 0(h 2 ) in H 1 , (^(h 1 ) in H 2 , and 0(h 2 ) in L 2 norms. 
The convergence rate for \\Q n (Vu ■ n e ) — u n \\t, is 0(h 2 ). Also, the same problem is 
tested for k = 1. The results are reported in Table 7.2. It indicates that the WG 
solution is convergent with rate 0(h 3 ) in if 1 , 0(h 2 ) in _ff 2 , and 0(h A ) in L 2 norms. 
We note that the L 2 error is convergent at order 4, two orders higher than that of 
k = 0, confirming the sharpness of Theorem 6.6. Moreover the convergence rate for 
||<9n(Vu • n e ) - u n \\b is OO 3 ), for k = 1. 



Table 7.2 

Example 1. Convergence rate element P%(T) — P2(e) (fc = 1) 



/l 


||u - 1*0 1 


{ \\u h - Qhu\\\ 


||u - Qou\\ 


||Q„(Vm • n e ) - u n \\ b 


1/4 


1.5888c-03 


1.5888e-02 


1.5751e-04 


1.7898c-03 


1/8 


2.6787c-04 


4.7921c-03 


1.3887e-05 


2.6200c-04 


1/16 


3.8354c-05 


1.2963c-03 


1.0006e-06 


3.4742c-05 


1/32 


5.0893c-06 


3.3568c-04 


6.6590c-08 


4.4563c-06 


1/64 


6.5373e-07 


8.5314c-05 


4.2842e-09 


5.6344e-07 


1/128 


8.2783c-08 


2.1499c-05 


2.7341e-10 


7.0798c-08 


Conv.Rate 


2.8597 


1.9152 


3.8450 


2.9336 



7.2. Example 2. In this problem, we set f2 = (0, l) 2 and the exact solution: 

u = sin(7rx) sin(7ry), 

with 

„ du . „ 
u an = 0, — ^ 0. 
on 

Boundary conditions and / are given according to the equation (7.1)-(7.3). 

Again, the uniform triangular mesh is used in the experiment. Table 7.3 shows 
that the convergence rate for C°-WG solutions in H 1 , H 2 and L 2 norms is 0(h 2 ), 
O(h) and 0(h 2 ), respectively. 



Table 7.3 

Example 2. Convergence rate for element P2(T) — Pi(e) (k = 0). 



h 


\\u - u \\i 


\\\u h - Qhu\\\ 


||«o - Qou\\ 


||<9„(Vu • n e ) - u n \\b 


1/4 


6.1653c-01 


5.5381 


1.2978e-01 


2.7515c-01 


1/8 


1.4737e-01 


2.7431 


3.2219e-02 


6.8563e-02 


1/16 


3.6122c-02 


1.3640 


7.9854c-03 


1.6489e-02 


1/32 


8.9758c-03 


6.8082c-01 


1.9899e-03 


4.0589e-03 


1/64 


2.2403c-03 


3.4024c-01 


4.9703e-04 


1.0102c-03 


1/128 


5.5983c-04 


1.7010c-01 


1.2423c-04 


2.5224e-04 


Conv.Rate 


2.0186 


1.0046 


2.0058 


2.0209 



7.3. Example 3. The exact solution is chosen as 

u = sin(7rx) cos(7ry), 
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with nonhomogeneous boundary conditions. 

Table 7.4 shows that the convergence rate for C°-WG solutions in H 1 , H 2 and 
L 2 norms is 0(h 2 ), 0(h), and 0(h 2 ), respectively. 



Table 7.4 

Example 3. Convergence rate for element Pi{T) — Pi(e) (k = 0). 



h 


\\u - u \\i 


\\\u h - Qhu\\\ 


||m - Qou\\ 


\\Q n (Vu ■ n e ) - u n \\ b 


1/4 


2.7134c-01 


4.3389 


2.8817e-02 


5.9389c-01 


1/8 


5.6175c-02 


2.4888 


5.8917e-03 


2.0490e-01 


1/16 


1.3236e-02 


1.3196 


1.3285e-03 


5.9347e-02 


1/32 


3.2856c-03 


6.7374c-01 


3.2089e-04 


1.5585c-02 


1/64 


8.2159e-04 


3.3917c-01 


7.9441c-05 


3.9554c-03 


1/128 


2.0553c-04 


1.6994e-01 


1.9812e-05 


9.9329e-04 


Conv.Rate 


2.0608 


9.4191c-01 


2.0916 


1.8609 



7.4. Example 4. In the final example, we test the a case where the exact so- 
lution has a low regularity in the domain f2 = (0,1) 2 . The exact solution is given 

by 

o/o / . 3$ . 9 

u = r 1 sin 3 sin - 

V 2 2 

where (r, 6) are the polar coordinates. It is known that u 6 H 2 b (£l). The perfor- 
mance for C° weak Galerkin finite element approximations for element PiiT) — P\(e) 
(k = 0) is reported in Table 7.5. The convergence rates in f/^-norm, H 2 — norm, and 
edge-based L 2 -norm are seen as 0(h 1A ) 1 0(h aA7 ) 7 and 0(h 1A ). The corresponding 
theoretical prediction has the order of 0(h 15 ), O(h 5 ), and 0(h 15 ). We believe that 
the numerical results are in consistency with the theory. Table 7.5 indicates that the 
numerical convergence rate in the standard L 2 is of order 0(h 188 ), which exceeds the 
theoretical prediction of 0(h 15 ). 

Table 7.6 contains some numerical results for the weak Galerkin element Ps(T) — 
P2(e) (k = 1). The convergence rates in i? 1 -norm, H 2 — norm, and edge-based L 2 - 
norm are seen as 0(h 15 ), O(h 5 ) 7 and 0(h 15 ), which are completely in consistency 
with the theory. For the element-based L 2 error, Table 7.6 indicates a numerical con- 
vergence rate of order 0(h 2A9 ), which is also consistent with the theoretical prediction 
ofO^ 2 - 5 ). 

Appendix A. A mass-preserving Scott-Zhang operator. We will prove the 
existence of an interpolation Q used in (4.1) and in the previous section, which is a 
special Scott-Zhang operator [21]. The new Scott-Zhang operator preserves the mass 
on each element and on each face, of four orders and three orders less, respectively, 
when interpolating H 1 (Q) functions to the finite element Vh functions. We shall derive 
the optimal-order approximation properties for the interpolation in the section, which 
leads to a quasi-optimal convergence of the weak Galerkin finite element method (3.6). 

The original Scott-Zhang operator maps u G i/ 1 (17) functions to C°-Lagrange 
finite element functions, preserving the zero boundary condition if u € H 1 ^). It is 
an Lagrange interpolation. All the Lagrange nodes ([4]) on one element are classified 
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Table 7.5 

Example 4- Convergence rate for element Pi{T) — -Pi(e) (fc = 0). 



h 


\\u - U \\l 


\\\u h - Qhu\\\ 


||u - Qou\\ 


||<9„(Vu • n e ) - u n \\b 


1/4 


3.1965c-02 


9.0667e-01 


3.3386e-03 


1.5615e-01 


1/8 


1.3596c-02 


6.8589e-01 


1.1209e-03 


6.2562e-02 


1/16 


5.1368c-03 


4.9952c-01 


3.1392c-04 


2.3370c-02 


1/32 


1.8697c-03 


3.5808c-01 


8.2158e-05 


8.4733e-03 


1/64 


6.7020c-04 


2.5488c-01 


2.0925c-05 


3.0321c-03 


1/128 


2.3855e-04 


1.8081c-01 


5.2718c-06 


1.0784c-03 


Conv.Rate 


1.4233 


4.6844c-01 


1.8767 


1.4415 



Table 7.6 

Example 4- Convergence rate for element P3(T) — ft(e) (fc = 1). 



h 


||w-«o||i 


t \\\u h - Qhu\\\ 


||m - Qou\\ 


L ||Qn(Vw • n e ) - u n \\ b 


1/4 


2.5197c-02 


5.0303e-01 


1.3671c-03 


4.7712c-02 


1/8 


8.9650c-03 


3.5619e-01 


2.4629c-04 


1.6900e-02 


1/16 


3.1718c-03 


2.5190c-01 


4.3679e-05 


5.9764c-03 


1/32 


1.1215c-03 


1.7812c-01 


7.7825e-06 


2.1130c-03 


1/64 


3.9652e-04 


1.2595c-01 


1.3812e-06 


7.4708c-04 


1/128 


1.4019c-04 


8.9063e-02 


2.4431c-07 


2.6413c-04 


Conv.Rate 


1.4984 


4.9966e-01 


2.4907 


1.4995 



into three types: 

corner nodes Cj : 3 vertex nodes in 2D, or all edge nodes in 3D, 

middle nodes mj : all mid-edge nodes in 2D, or mid-triangle nodes in 3D, 

internal nodes ij : all internal nodes in the triangle/tetrahedra. 

The three types of nodes are illustrated in Figures A.l and A. 2. In simple words, {cj} 
are nodes shared by possibly more than two elements, {rrij} are nodes shared by no 
more than two elements, and {ij} are nodes internal to one element. 

A Lagrange nodal basis function <f>j is a Pk+2 polynomial which assumes value 1 
at one node Cj, but vanishes at all other dim P^ +2 — 1 nodes. For example, a P| nodal 
basis function on the reference triangle {0 < x, y, 1 — x — y < 1}, at node (1/4, 0), c.f. 
Figure A.3, is 

, An A( x x(l-x-y)(3/4-x-y)(2/A-x-y) 
(A.l) <j> 2 {x,y)- 



(1/4)(1 - 1/4 - 0)(3/4 - 1/4 - 0)(2/4 - 1/4 - 0) ' 

The restriction of a nodal basis <j>j on a lower dimensional simplex, a triangle or an 
edge or a vertex, is also a nodal basis function on that lower dimensional finite element. 
For example, this node basis function (A.l) is the restriction of the following 3D nodal 
basis function (at node (1/4, 0, 0) on tetrahedron {0 < x, y, z, 1 — x — y — z < 1}) on 
the reference triangle, 

(A2] 6 (x v z)= x(l-x-y- z)(3/4 -x-y- z)(2/4 - x - y - z) 
y ~> ^ ' y ' j (1/4)(1- 1/4-0- 0)(3/4- 1/4-0 -0)(2/4- 1/4- 0-0)" 
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Fig. A.l. Averaging patches (Cj) in 2D (an edge or a triangle), for corner nodes (cj), middle 
nodes (mj) and internal nodes (ij), in 2D. 










\ • V. 












Fig. A. 2. Averaging patches (Cj) in 3D (a triangle or a tetrahedron), for corner nodes (cj), 
middle nodes (mj) and internal nodes (ij), in 3D. 



The restriction of 2D basis function </> 2 in (A.l) in ID is, c.f. Figure A. 3, 

x(l -x)(3/4-x)(2/4-x) 



(A.3) 



Pj' 



(x) 



(l/4)(l-l/4)(3/4-l/4)(2/4-l/4)- 





d)j> (in ID) 
• ® • • • 

Fig. A. 3. A 2D nodal basis <f>2 (A.l) and its restriction in ID, <f>ji (A.3). 

On each element T (an edge, a triangle, or a tetrahedron), the Pk Lagrange basis 
{4>j} has a dual basis {tpj € P^}, satisfying 



(A.4) 



In other words, if writing {tpj} as linear combinations of Lagrange basis {4>j}, the 
coefficients are simply the inverse matrix of the mass matrix, the L 2 -matrix of {<j>j}. 



1G 



For example, the dual basis function ip 2 for the nodal basis function <j> 2 in (A.l) (2D) 
is 



V4 2DI (x, y) 



2835 



12285 



8505 



x 3 + 8505a; 2 ?/ + 



8505 



4 4 2 

- 1890a; 4 - 5670a; 3 ?/ - 5670a; 2 ?/ 2 



-xy 



1890a;y d 



We can compute the dual of ipj> in (A. 3) in ID to get 

.rid, , 485 2865 21105 2 13615 
(A.5) 4 ](*) = -— + —*-—*» + . 



32 



11655 

IT' 



485 , 64225 , 345 

4>ci H <Pi H <; 

128 16384^ 1024 s 



4255 
16384' 



85 
128* 



where <pi is the nodal basis on [0,1] at Xi = i/A, i = 0,1,2,3,4. The dual function 
ip^ D \x) in (A.5) is plotted in Figure A. 4. Similarly we can compute the dual basis 
function for (A. 2) in 3D. We note that both Lagrange nodal basis and its dual basis 
are affine invariant. That is, the Lagrange basis on the reference triangle is also the 
Lagrange basis on a general triangle after an affine mapping. For simplicity, we use 
the same notations <j>j and ipf D ^ for the nodal basis and the dual basis functions on 
the reference triangle and on a general triangle. 




Fig. A. 4. A ID nodal basis <j)j> (A. 3) and its dual basis in ID, if>^, D ^ (A.5), 4>j'i'^ 1 / D ^ = 1- 



We now define the Scott-Zhang interpolation operator: 

Qo : tf^ft) -> V},, 

where Vh is the C°-Pk+2 finite element space defined in (3.2). Q v is defined by the 
nodal values at three types nodes. 

1. For each corner node Cj (shared by possibly more than two elements), we 
select one boundary (d— l)-dimcnsional simplex Cj if Cj is a boundary node, 
or any one (d — l)-dimensional face simplex Cj on which the node is, as c^-'s 
averaging patch. C.f. Figure A.l, the boundary node Cj has a boundary edge 
Cj, while the corner node Cj" can choose any one of four edges passing it, as 
its averaging patch. In Figure A. 2, a corner node Cj has a triangle Cj as its 
averaging patch. In both 2D and 3D, we use a same definition 

(A.6) Q v( Cj )= I S d - 1)D] v(x)dx. 

Jc 3 



2. For each middle node m,j, the averaging patch is the unique (d— l)-dimcnsional 
simplex Cj containing m,j, see Figures A.l and A. 2. The interpolated nodal 
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value is then determined by the unique solution of linear equations: 
(A.7) 

/ ( Qov(mji)4>j>(x) + ^2 Qov(cj)4>j(x) - v(x)^Jpi(x)dx = 

for all degree (k + 2 — d) polynomials pi on (d — l)-simplex Cj , where Qov(cj ) 
is defined in (A. 6). In 2D, c.f. Figure A.l, after we determine the nodal 
values at the two end points (cj), we determine the middle-edge points' nodal 
value by (A.7). 

3. After determine all nodal values on the surface of each element, we define the 
interpolation inside the element: The nodal values at internal nodes (Qo{ij")) 
are determined by the unique solution of the following linear equations 



(A.8) / ( X! Qov(ij")(pj"^Pidx 



= (v - Qov(mj>)<j)j> - ^2 Qov{cj)4>j^pidx 

for all degree (k + 1 — d) polynomials p, on d-simplex Cj. 
By (A.6)-(A.8), the (refined) Scott-Zhang interpolation is 

(A.9) Qqv = ^2 Qov(xj)(f>j(x), 

where Mh is the set of all C°-Pk+2 Lagrange nodes of triangulation Th ■ 

Remark A.l. If all corner nodes have selected a same averaging patch Cj as 
the unique patch for the middle nodes on the patch, then the solution of (A.7) is the 
L 2 -projection, i.e., 

(A.10) Qov(mj>) = [ i>f, d ~ 1)D] v{x)dx. 

JCj 

In fact, (A.10) is the definition of the original Scott-Zhang operator in [21]. In the 
same fashion, if all patches of the face nodes are face (d — l)-simplexes of Cj, then 
the internal nodal values are exactly that of the L 2 -projection on Cj, i.e., the solution 
of (A.8) satisfies 

(A.11) Qow(V') = / ^ D] v{x)dx. 

Jc 3 

But if there are more than one triangle or tetrahedron in Th, some Cj must be from 
neighboring elements. So (A.10) and (A.11) can not be satisfied in general. Otherwise 
the Scott-Zhang operator would preserve mass of order (fc + 2) both on an element and 
on its faces. 

Lemma A.l. The Scott-Zhang interpolation operator (A.9) is well-defined, i.e., 
the linear systems of equations (A.7) and (A.8) both have unique solutions. 

Proof. For the linear system of equations (A.7), we change the P^+l-d basis 
functions {{pi = 1, x, x h } when d = 2, or {pi = 1, x, y, x 2 , xy, y h ~ 1 } when d = 3) 



18 



uniquely as linear combinations of the Lagrange basis functions on the subinterval Cj 
(d = 2) or the subtriangle Cj (d = 3), c.f. Figure A.5. 



(A.12) 



The nodal basis functions on a simplex Cj and its subsimplcx Cj differ by a bubble 
functions: 



(A.13) 



6(x) 



where 6(x) is the bubble functions assuming on the boundary of Cj. For example, 
c.f. Figure A. 5, when d — 2 and Cj = [0, 1], 

b(x) = x(l — x), 

(ar-2/4)(a;-3/4) 



6 2 (a;) = <^(x 



(l/4-2/4)(l/4-3/4)' 
6(a0 



6(1/4)' 



d= 2: 



d = 3: 



(z-2/4)(a;-3/4) 
(l/4-2/4)(l/4-3/4) 



2 = 4>\b{x)/b(l/A) 



C? 



N 



Fig. A. 5. Lagrange nodal basis <j> a - on subinterval (d = 2) or subtriangle (d = 3), c.f., (A.13). 

By the change of basis, (A.12) and (A.13), the linear system (A. 7) is equivalent 
to the following weighted-mass linear systems: 

(A.14) Qov(m r ) ^ r ((> s i dx= vftdx - ^ Q v(cj) / fafidx, 

i=l,2,...,dim(P fc rf - 2 1 _ tI ). 

The coefficient matrix in (A.14) is the mass matrix on the subsimplex Cj (c.f. Figure 
A. 5) with a positive weight: 



a i,j' = / <Pi<Pj'd x = / <t>t4>jw(x)dx, 



where 



w(x) = 



*>(*) 
6(xj-/) 



> in interior(Cj). 



As the Lagrange basis {0| } (on the subsimplex) are linearly independent, the mass 
(with weight) matrix in (A.14) is invertible, and the equivalent linear system (A. 7) 
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has a unique solution too. For example, when d = 2 and k = 2, the coefficient matrix 
in (A. 14) and its inverse are 

/152/315 -16/63 8/63 \ _1 / 2655/1024 75/64 -225/1024\ 

-4/21 18/35 -4/21 = 225/256 45/16 225/256 
\ 8/63 -16/63 152/315/ \-225/1024 75/64 2655/1024/ 

By the same argument, lifting the space dimension by 1, we can show that (A. 8) has 
a unique solution too. In fact, the system (A. 8) when d — 2 is the same system (A. 7) 
with d = 3 there. □ 

Lemma A. 2. If d=2, the Scott-Zhang interpolation operator (A. 9) preserves the 
volume mass of order k — 1 and the face mass of order k, i.e., 

(A.15) J (v- Qov) Pl dx = VT e Th, Pi <E PLi(T), 

(A.16) f (v - Q v) Pt dx = VT e £ h , Pi G P k {E), 

JE 

where £h is the set of edges in triangulation Th ■ 
Proof. By the construction (A. 7), we have 

(v-Q v)pidx = / (v- ^2 Qov{mj>)(t>j> - ^ Q v(cj)(f>j)pidx = 0. 

E E m jt eE c 3 GB 

That is (A.16). Here, because we have two missing dof (degrees of freedom) at the 
two end points of each edge, the polynomial degree in mass preservation is reduced by 
two, from (k + 2) to k. Similarly, (A.15) follows (A. 8). Here, the polynomial degree 
reduction is three as each triangle has three edges where the interpolation values are 
not free (not determined by (A. 8)). □ 

Lemma A. 3. If d = 3, the Scott-Zhang interpolation operator (A. 9) preserves the 
volume mass of order k — 2 and the face mass of order k — 1, i.e., 

(A.17) J (v - Q v) Pi dx - VT g T h , Pi G P fe 3 _ 2 (T), 

(A.18) f (v - Q v) Pi dx - VT G £ h , Pl G P^E), 

JE 

where £h is the set of face triangles in the tetrahedral grid Th- 

Proof. As each triangle E has three edges, where the interpolation is not deter- 
mined possibly by function value on neighboring triangles, we lose dof's on the three 
edges in the interpolation. That is, we lose three orders in face-mass conservation in 
(A. 7). (A.18) is simply another expression of (A. 7), as in the proof of Lemma A. 2. 
By (A. 8), (A.17) follows. Here the polynomial-degree deduction in mass conservation 
is 4, due to 4 face-triangles each tetrahedron. □ 

Remark A. 2. By Lemmas A. 2 and A. 3, the mass preservation on element and 
on faces is one order higher than the requirement (4.1), when d = 2. When d = 3, 
(A.17) and (A.18) imply (4.1). 
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Theorem A. 4. The Scott-Zhang interpolation operator (A. 9) is of optimal order 
in approximation, i.e., when k > — 1, 

(A.19) \\v-Qov\\+h\\v-Q v\\i < Ch k+3 \\v\\ k+3 V« G H k+3 (n). 

Further, when k > 0, 

(A.20) ( J2 h2 \ v - Qov\ 2 H 2 (T) ) 1/2 < C^ +3 |M| fe+3 V« G nH k+3 (Q). 

Ten 

Proof. The Scott-Zhang operator preserves a degree (k + 2) polynomial on the 
star union of an element T, 

St = Uy7 nT ^ T', T,T' G Th- 

That is, 

(A.21) Q v = v V«eP^(S T ), 

when k > —1. (A.21) is shown in three steps. First, by (A. 6), when v G Pk+2, the 
dual basis defines 

(A.22) Qov(cj) = v{cj), 

at all corner nodes. In the second step, by (A.22) and (A. 14), (A. 7) holds for pi = i[>* 
too where m,' are all middle nodes on Cj. That is, 

(A.23) f (v-J2 Qo(cj)4>j - J2 Qov(mj')4j')<l>idx = 0. 

By (A.22), as v G P k+2 , 

v - ^ Q (cj)<f>j = w c 6(x) for some w c G P^-d' 

where 6(x) is a bubble function, cf. (A. 13). For the middle node basis functions, we 
have also, cf. (A. 13), 

<j>y = (f>j„b(x)/b(xj>) Vrrij' G Cj. 
Thus (A.23) implies, where J2 m -, Qo v ( m j') ( t>j' — v m b(x) for some v m G Pk+2-d' 

/ («h - -y m )^ s 6(x)dx = 0, 

=> «b — w m = and w = Q W on Cj ■ 
Therefore, at all middle nodes, 
(A. 24) Q v(rrij') = v(m 3 <). 

In the third step, by (A. 8) and the same argument in the second step, 

Qo(V') = v {if')i 
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at all internal nodes. Thus Q v = t£ Pk+2- 

We then use the standard scaling argument (on the dual basis functions) and the 
Sobolev inequality, as in Theorem 3.1 of [21], it follows that 

\Qov\m( T ) < C\\v\\ m{ST) VveH'in). 

The above stability result leads directly to the optimal-order approximation (A. 19), 
following the standard argument (i.e., by (A. 21) and the existence of local Taylor 
polynomials, c.f. for example, [4]), as shown in Theorem 4.1 of [21]. We note again 
that the Scott-Zhang operator here is a refined version of the Scott-Zhang operator 
in [21]. After showing the local preservation of Pk+2 polynomials above, the proof of 
the theorem is the same as that in [21]. (A. 20) is (4.4) in [21], with p = q = 2, m = 2, 
and I — k + 3 there. 
□ 
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